jacobian <- function(x, prods, T, nn, v, Y, K, x_char, prodsMarket, numProdsTotal, marketStarts, marketEnds, oo, sharesum, marketForProducts, nD) { 

x<-as.matrix(x)

theta1 <- x[1:(K+1), 1]     
theta2 <- x[((K+1)+1):((nD+1)*(K+1)),]    
delta = x[((nD+1)*(K+1)+1):((nD+1)*(K+1)+numProdsTotal), 1]               
g = x[((nD+1)*(K+1)+numProdsTotal+1):nrow(x), 1]  
g <- as.matrix(g) 
  

expmu = matrix(0,nrow(x_char),ncol(v))
for (t in 1:T) { 
expmu[marketStarts[t]:marketEnds[t],] = exp(x_char[marketStarts[t]:marketEnds[t],]%*%diag(theta2[1:(K+1)])%*%v[1:(K+1),,t])
for (dd in 2:nD) {
expmu[marketStarts[t]:marketEnds[t],] = expmu[marketStarts[t]:marketEnds[t],] * exp(x_char[marketStarts[t]:marketEnds[t],]%*%diag(theta2[((K+1)+1):(2*(K+1))])%*%v[((dd-1)*(K+1)+1):(dd*(K+1)),,t])
}
}


expmeanval = as.matrix(exp(delta))
simEstShare<-ind_shnormMPEC(expmeanval,expmu, oo, sharesum, marketForProducts)
simShare<-simEstShare[,1:(ncol(simEstShare)-1)]
EstShare<-simEstShare[,ncol(simEstShare)]
EstShare<-as.matrix(EstShare)


ooo<-matrix(1,1,K+1)
  
# Evaluate the Gradients
dSdtheta2<-matrix(0,numProdsTotal,(nD*(K+1)))
dSddeltaDIAG<-matrix(0,numProdsTotal,prods)
dSddelta<-matrix(0,numProdsTotal,numProdsTotal)
  

for (t in 1:T) {
  index<-marketStarts[t]:marketEnds[t]
  ooo1<-matrix(1,prodsMarket[t],1)
  index_start <- marketStarts[t]
  index_end <- marketEnds[t]
  
  dSddeltaDIAG_old <- matrix(0, length(index),prodsMarket[t])
  dSddeltaDIAG_temp = dc_1_cpp (simShare, index_start, index_end, dSddeltaDIAG_old, nn )
  dSddeltaDIAG[index,1:prodsMarket[t]] = dSddeltaDIAG_temp
  dSddelta[index,index] = dSddeltaDIAG[index,1:prodsMarket[t]]
  
  dSdtheta2 = dc_2_cpp_new(dSdtheta2, simShare, ooo, ooo1, v[,,t], x_char, nn, index_start, index_end, nD, K, t) 
  
}


Ddelta = -solve(dSddelta)%*%dSdtheta2


return(Ddelta)
  
}  